fix bMax quality issues via objective offset / constant extraction - #4

1. Summary

$$ D_1=\int_{0}^{b_{\max}}\frac{b^{N+1}}{\prod_k \sqrt{\sigma_k + 2b^2}} db $$

$$ D_2=\int_{0}^{b_{\max}} b\prod_k \sqrt{\frac{2b^2}{\sigma_k^2+2b^2}} \sum_i w_i \exp \left( -\frac12 \sum_k \frac{x_{ik}^2} {\sigma_k^2+2b^2} \right) db $$

$$ D_3=\frac18 \sum_{i\ne j} w_i w_j \left[ 4b_{\max}^2 +\left(\gamma-1-2\ln 2b_{\max}\right)d_{ij}^2 +d_{ij}^2\ln d_{ij}^2 \right] $$

The problem is that $D_2$ and $D_3$ each carry a large constant of their own and those were not being removed.

2. Where the constants lives

2.1 The $D_2$ term

From calculateD2 / calculateP2 (gm_to_dirac_short.tpp):

$$ D_2 = -2\int_0^{b_{\max}} \rho(b) \sum_{i=1}^{L} w_i e^{z_i(b)} db $$

with the width-dependent prefactor and exponent

$$ \rho(b) = b\prod_{k=1}^{N}\sqrt{\frac{2b^2}{\sigma_k^2+2b^2}}, \qquad z_i(b) = -\frac{1}{2}\sum_{k=1}^{N}\frac{x_{ik}^2}{\sigma_k^2+2b^2}. $$

The position dependence enters only through $z_i(b)$. Now observe the limit: as $b\to\infty$ the denominators $\sigma_k^2+2b^2$ blow up, so $z_i(b)\to 0$ and $e^{z_i(b)}\to 1$ for every point, regardless of where it is. Meanwhile $\rho(b)\to b$, which grows without bound.

So the integrand tends to $b\sum_i w_i = bW$: a linearly growing quantity carrying zero information about the configuration. Integrating it over $[0,b_{\max}]$ contributes $\approx -b_{\max}^2 W$ of pure dead weight.

Making the split explicit with the exact identity $e^{z} = 1 + expm1(z)$ (where $expm1(z) = e^{z} - 1$):

$$ \sum_{i} w_i e^{z_i(b)} = \underbrace{\sum_i w_i}_{\textstyle W} + \sum_i w_i expm1{\big(z_i(b)\big)} $$

$$ D_2 = \underbrace{-2W \int_0^{b_{\max}} \rho(b) db}_{\text{constant in }x} \quad - \quad \underbrace{2\int_0^{b_{\max}}\rho(b)\sum_i w_i \; {expm1}{\big(z_i(b)\big)} db}_{\text{depends on }x} $$

The constant piece is already computable from existing code. Writing $\rho$ in terms of calculateP1:

$$ \rho(b) \;=\; b\prod_k \frac{\sqrt{2},b}{\sqrt{\sigma_k^2+2b^2}} \;=\; 2^{N/2}\underbrace{\,b\prod_k \frac{b}{\sqrt{\sigma_k^2+2b^2}}\,}_{\textstyle P_1(b)} \;\;\Longrightarrow\;\; \int_0^{b_{\max}}\!\!\rho\,db \;=\; 2^{N/2} D_1 . $$

So the $D_2$ constant is exactly $-2W\,2^{N/2}D_1$, where D1 and twoPiNHalf $=2^{N/2}$ are already members of GMToDiracBaseOptimizationParams.

Note: the member is named twoPiNHalf but is initialised to $2^{N/2}$, not $(2\pi)^{N/2}$. The $(2\pi)^{N/2}$ line is commented out in the source (gm_to_dirac_optimization_params.h, line 80-81) . The name is misleading, the value is what this derivation uses.

2.2 The $D_3$ term

From calculateD3, the pairwise self term is

$$ D_3=\frac18 \sum_{i\ne j} w_i w_j \left[ 4b_{\max}^2 +\left(\gamma-1-2\ln 2b_{\max}\right)d_{ij}^2 +d_{ij}^2\ln d_{ij}^2 \right], \qquad d_{ij}^2=\lVert x_i-x_j\rVert^2 $$

The leading $4b_{\max}^2$ has no $d_{ij}$ in it at all. Summed over every ordered pair:

$$ \sum_{i,j}\frac{w_iw_j}{8}\cdot 4b_{\max}^2 \;=\; \frac{b_{\max}^2}{2}\Big(\sum_i w_i\Big)\Big(\sum_j w_j\Big) \;=\; \frac{1}{2}b_{\max}^2W^2 . $$

2.3 The total offset, in closed form

$$ \boxed{\;C(b_{\max}) \;=\; -2W\,2^{N/2}D_1(b_{\max}) \;+\; \tfrac{1}{2}\,b_{\max}^2W^2\;} $$

The optimizer now optimises $D_2 + D_3$ without the constants. The public distance function still returns the old value.

3. Results

3.1 Standard Normal Samples

3.2 Correlated Samples